Clustering of inertial particles in a turbulent flow 
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We analyzed formation of small-scale inhomogeneities of particle spatial distribution (particle 
clustering) in a turbulent flow. The particle clustering is a consequence of a spontaneous breakdown 
of their homogeneous space distribution, and is caused by a combined effect of the particle inertia 
and a finite correlation time of the turbulent velocity field. Theory of the particle clustering is 
extended to the case when the particle Stokes time is larger than the Kolmogorov time scale, but 
is much smaller than the correlation time at the integral scale of turbulence. The criterion of the 
clustering instability is obtained. Applications of the analyzed effects to the dynamics of inertial 
particles in industrial turbulent flows are discussed. 

PACS numbers: 47.27.Qb 
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I. INTRODUCTION 

It is generally believed that turbulencepromotes mix- 
ing (see, e.g., 0,111110113). However, 
laboratory experiments and observations in the atmo- 
spheric turbulence show formation of long-living inho- 
mogeneities in concentration distribution of small iner- 
tial particles and droplets in turbulent fluid flows (see, 
e -g-; Uil EH E3 ) ■ The origin of these inhomogeneities is 
not always clear but there influence on the mixing can 
be hardly overestimated. 

The goal of this study is to analyze the particle-air 
interaction leading to the formation of strong inhomo- 
geneities of particle distribution, referred to as particle 
clustering (see |131 ll 1|. and references therein). Particle 
clustering is a consequence of a spontaneous breakdown 
of their homogeneous space distribution. As a result of 
the nonlinear stage of clustering, the local density of par- 
ticles may rise by orders of magnitude and strongly in- 
crease the probability of parti cle-p article collisions. 

It was suggested in IIH fl6| that the main reason 
for the particle clustering is their inertia: the particles 
inside the turbulent eddies are carried out to the bound- 
ary regions between the eddies by the inertial forces. 
This mechanism of the preferential concentration acts in 
all scales of turbulence, increasing toward small scales. 
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Later, this was contested in Refs. 0, E3 using the so- 
called "Kraichnan model" 0] °f turbulent advection 
by the delta-correlated in time random velocity field, 
whereby the clustering instability did not occur. How- 
ever, it was shown in Ref. 0[ that accounting for a finite- 
ness of correlation time of the fluid velocity field results 
in the clustering instability of heavy particles. The the- 
ory of the clustering instability of the inertial particles 
advected by a turbulent velocity field was developed in 
Ref. 0| and applied to the dynamics of aerosols in the 
turbulent atmosphere. 

In the present study we extended the theory of parti- 
cle clustering to the case when the particle Stokes time 
is larger than the Kolmogorov time scale, but is much 
smaller than the correlation time at the integral scale 
of turbulence. The paper is organized as follows. In 
Sec. [n] we present a qualitative analysis of the clustering 
instability that causes formation of particle clusters in a 
turbulent flow. In Sec. II I II we evaluate the characteris- 
tic parameters that affect clusters formation as a result 
of the clustering instability: the particle response time 
t p , the Kolmogorov micro-scale time r^, the characteris- 
tic velocity of small and large particles in turbulent fluid; 
the degree of compressibility of the particle velocity field, 
cr„(a), where a is the particle size. In Sec. lIIIBl we study 
the velocity of small and large particles in the turbulent 
fluid, required for evaluation of the effect of turbulent 
diffusion in these two regimes. In Sec. IIVI we present a 
quantitative analysis for the clustering instability of the 
second moment of particle number density. This allows 
us to generalize the criterion of the clustering instability, 
obtained in Ref. [TJ|]. Finally, in Sec. E] we overview the 
nonlinear effects which lead to saturation of the cluster- 
ing instability and determine the particle number density 
in the cluster. In Sec. we also present numerical esti- 
mates for droplet dynamics for the conditions pertinent 
to turbulence in diesel engines. 



2 



II. CLUSTERING OF PARTICLES IN 
TURBULENT GAS 

A. Basic Equations 

To analyze dynamics of particles we use the standard 
continuous media approximation, introducing the num- 
ber density n(t, r) of spherical particles with radius a. 

The particles are advected by a turbulent velocity field 
u(i,r). Since typically the velocity of the carrier gas is 
much smaller than the sound velocity, the incompressibil- 
ity constrain divu(/j, r) = is applicable. The particle 
material density p p is much larger than the density p of 
the ambient fluid. For hea vy p articles v(t, r) ^ u(t, r) 
due to the particle inertia |20i I2TL |22|. Therefore, the 
compressibility of the particle veloci ty fi eld v(i, r) must 
be taken into account (see 0, 0, 0, 0), since the 
growth rate of the clustering instability, 7, is proportional 
to (|div v(i, r)| 2 ), where (•) denotes ensemble average. 

Let 0(i, r) be the deviation of the particle number den- 
sity n(t, r) from its uniform mean value n = (n): 

G(t,r) = n(t,r)-fi. (2.1) 

The pair correlation function of 9(<, r) is defined as 

$(R, r, t) = (Q(t, r + R)9(t, r)> . (2.2) 

For the sake of simplicity we will consider only a spa- 
tially homogeneous, isotropic case when <&(R, r, t) de- 
pends only on the separation distance R and time t: 



$(t,R,r) -> <f>(t,R) 



(2.3) 



Denote the probability of the pair collisions between par- 
ticles as p(t) which can be expressed as 



p(t) oc c <M + 



(KM)] 2 ) 



}- 


L , $(f,o)l 







(2.4) 



Obviously, a large increase of &(t, R) above the level of 
n 2 leads to a strong grows in the frequency of the particle 
collisions. 

In the analytical treatment of the problem we use the 
standard equation for n(t,r): 



dn(t, r) 
dt 



+ V • [n(t, r)v(t, r)] = DAn(t, r), (2.5) 



where D is the coefficient of molecular (Brownian) diffu- 
sion. The equation for 0(i,r) follows from Eq. 1)2. 5fl : 



dO{t,r) 
dt 



+ [v(t,r)-V]9(<,r) (2.6) 
= -6(t,r)divv(t,r) + DAB(t,r) . 



Here we neglected the term oc ndivv, describing the ef- 
fect of an external source of fluctuations. It was shown in 
Ref. 0| that this effect is usually much smaller than the 
effect of self-excitation of fluctuations of particle number 
density. 



One can use Eq. i|2.6|) to derive equation for <&(/;, R) by 
averaging the equation for 6(i, r + R)0(i, r) over statis- 
tics of the advected turbulent velocity field v(rj, r). In 
general this procedure is quite involved even for simple 
models of the advecting velocity fields, see, e.g., Ref. [l4| . 
Nevertheless, the qualitative understanding of the un- 
derlying physics of the clustering instability, leading to 
both, the exponential growth of R) and its nonlin- 
ear saturation, can be elucidated by a more simple and 
transparent analysis, that is presented below. 



B. Qualitative analysis of the clustering instability 

1. On Richardson-Kolmogorov cascade theory of turbulence 

In our discussion wo will use the well known 
Richardson-Kolmogorov cascade theory of turbulence 
(see, e.g., Refs. Q, 123,0). For the large Reynolds num- 
bers Tie ^> 1 the characteristic scale L of energy injection 
(outer scale) is much larger than the length of the dissi- 
pation scales (viscous scale rj) L rj. In the so-called 
inertial interval of scales, where L > r > 77, the statistics 
of turbulence within the Kolmogorov theory is governed 
by the only dimensional parameter, e, the rate of the tur- 
bulent energy dissipation. Then, the velocity u(r) and 
the energy of turbulent motion E[r) at the characteristic 
scale r (referred below as r-eddies) may be found by the 
dimensional reasoning: 

u(r) « [erf\ E(r) = \p[u{r)f . (2.7) 

Similarly, the turnover time of r-eddies, r(r), which is of 
the order of their life time, may be estimated as 



r(r) ~ r/u(r) w e ' r 



1/3J2/3 



(2.8) 



To elucidate the clustering instability let us consider a 
cluster of particles with a characteristic scale t moving 
with the velocity V c i(t). The scale £ is a parameter which 
governs the growth rate of the clustering instability, 7. 
It sets the bounds for two distinct intervals of scales: 
L > r > £ and £ > r > rj. However, if the size of particles 
a is larger than the viscous scale 77, the second range of 
scales becomes £ > r > a, because we cannot consider 
scales which are smaller than the size of particles. Large 
r-eddies with r > £ sweep the ^-cluster as a whole and de- 
termine the value of V c i(i). This results in the diffusion 
of the clusters, and eventually affects their distribution 
in a turbulent flow. The small r-eddies determine the dy- 
namics of particles inside the cluster. The role of these 
eddies is multifold. First, they lead to the turbulent dif- 
fusion of the particles within the scale of a cluster size. 
Second, due to the particle inertia they tend to accumu- 
late particles in the regions with small vorticity, which 
leads to the preferential concentration of the particles. 
Third, the particle inertia also causes a transport of fluc- 
tuations of particle number density from smaller scales 
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to larger scales, i.e., in regions with larger turbulent dif- 
fusion. The latter can decrease the growth rate of the 
clustering instability (see below). Thus, the clustering is 
determined by the competition between these three pro- 
cesses. 



2. Effect of turbulent diffusion 

We consider Eq. (|2.6I) for 0(t, r) in the coordinate sys- 
tem co-moving with the ^-cluster and assume v(t, r) w 
u(i,r). In this reference frame the advectived velocity 
v(i,r) should be replaced by [v(i, r) — V c i(t)]. In par- 
ticular, the advection term in Eq. I|2.6[l takes the form 



Adv={[v(t,r)-V c i(t)]-V}0(t,r) 



(2.9) 



Averaging this term over statistics of turbulent velocity 
field leads to the turbulent diffusion of particles. It is well 
known (see, e.g. Ref. 0) that the turbulent diffusion can 
be modelled by the renormalization of the diffusion term 
in the right hand side (RHS) of Eq. (|2~^|l : 



D=> D + D, 



(2.10) 



where D T is the turbulent diffusion coefficient. The main 
contribution to the advected velocity in Eq. (|2.9J) is due 
to the velocity of ^-eddies, v(£). Therefore D T becomes 
a function of scale £, D T D T (£), and can be estimated 
from the parameters of ^-eddies, using dimensional rea- 
soning: 



D T (£)*~£v(£) 



(2.11) 



Here we used the commonly accepted [23|, |24| numeri- 
cal factor i for the turbulent diffusion coefficient in an 
isotropic turbulence. Now the part of Eq. i|2.6[l describing 
the turbulent diffusion may be written as 



dG(t,r) 

at 



-£> T (£)Ae(i,r) 



(2.12) 



During the linear stage of the cluster evolution the 
particle distribution inside the cluster does not change. 
Therefore, the function 0(i, r) in the expression for the 
correlation function $ can be factorized: 



3. Effect of particles inertia 

For the qualitative analysis of the particle inertia we 
consider Eq. I|2.6|) . written in the co-moving reference 
frame, taking into account only inertia term: 



dQ(t,r) 
dt 



= -0(t,r)div[v(t,r)-V cl (t)] 



(2.15) 



As before we can factorize the function Q(t, r) according 
to Eq. I|2.13[) and simplify the partial differential equa- 
tion (|2.15l) . reducing it to the differential equation for 
the cluster amplitude Ai(t): 



dA e (t) 
dt 



-A e (t)b e (t) 



(2.16) 



Here we neglected the r-dependence of the divergence 
term inside the cluster in the RHS of Eq. (|2.15[) . so that 
this term becomes a function of t only: 



div[v(i,r)-V c i(i)]=>6*(t) . 



(2.17) 



In the turbulent flow field the function bg(t) may be 
considered as a random process with some correlation 
time Tb, which will be evaluated below. Since the insta- 
bility of ^-clusters is caused by ^-eddies, the correlation 
time Tb can be estimated as the turnover time of the par- 
ticle velocity field of ^-eddies: 



n » £/v{£) 



(2.18) 



The evaluation of the mean square value of be(t) re- 
quires a more careful consideration. One cannot estimate 
divv(r) by the dimensional reasoning as v(r)/r. Indeed, 
in the incompressible flow divv(r) = 0. To elucidate this 
issue we introduce a dimensionless parameter a v , a de- 
gree of compressibility of the velocity field of particles in 
f-clusters, V£(t, r), defined by 



<[divv,] 2 )/(|V x v,| 2 ) 



(2.19) 



This parameter may be of the order of 1 (see Refs. 0, 
At the moment we can evaluate be(t) via this yet 
unknown parameter a v as follows: 



(&/(*)> =0, (bj(t))^<v(£)/£f 



(2.20) 



Let us show that the stochastic differential equation 
(|2.16l) results in an exponential growth of (Af(t)), i.e. 
in the instability. Indeed, the solution of this equation 
reads: 



( 2 



(2.13) 



®{t,r)^A t (t)6 t {r) 
Then, Eq. TFFfy yields 

A}{t) = A 2 (0)exp[-7dift] , 7dif = 

i - 

(2.14) 

where the Laplace operator is evaluated as l/£ 2 and we 
used Eq. lE~nj). 



2D T (£) 2v(£) 



A e (t) = A e (0) exp[-/(i)] , I(t) = J b e (T)dr . (2.21) 

'o 

Integral I(t) may be written as the sum of integrals I n 
over small time intervals 7~(,: 



I(t) 



N 



(n+l)r b 

7„ = J be(r)dT, 

UTb 



N = — . (2.22) 

n 



4 



By definition, T& is the correlation time of the random 
process b^(t). Therefore the integrals I n can be consid- 
ered as independent random variables. According to the 
central limit theorem, the sum of a large number of sta- 
tistically independent random variables is distributed as 
the Gaussian random variable. Therefore the total inte- 
gral I(t) can be estimated as 



(2.23) 



Here £ is a Gaussian random variable with zero mean 
and unit variance. The probability density function of £ 
reads 



v(C) 



1 



'2it 



exp 




(2.24) 



Using Eqs. {223}, (223) and lt2~2~2"|) we obtain 



(Aj(t)) ~ A, 2 (0) / exp[-2^/N(lZ)t]V(0dC 



where the growth rate 7j„ is given by 



7in 



n 



2 n (b 2 e (t)) 



2a v v{€) 
I 



(2.25) 



(2.26) 



In the last estimate we used Eqs. H2.18|l and 12.20JI . It is 
clearly seen that the source of the instability is a nonzero 
value of a v , i.e. a compressibility of the particle velocity 
field, v(t, r). 

In the case of the small enough particles, t p < r(r]), 
where r(rj) is the turnover time of the smallest, Kol- 
mogorov micro-scale eddies. In this case all the particles 
are almost fully involved in turbulent motion, and one 
concludes that u(t, r) sa v(t, r) and v(£) ~ u{t). Also 
we will see below that in this case the largest value of 
o~ v is attained for I w r\. Hence for r p < r(r\) the most 
unstable clusters are those with I « 77. In summary: 

• In the case r p < r(rf), the characteristic scale of the 
most unstable clusters of small enough particles is of the 
order of Kolmogorov micro-scale of turbulence, n. 

• The characteristic growth rate of the clustering insta- 
bility is of the order of the turnover frequency of 77-eddies, 
1/t( V ). 

• The particle clusters are unstable for heavy enough 
particles, such that the degree of compressibility a v (rj) 
of their effective advective velocit y fi eld exceeds some 
threshold value cr cr « 0.3 (see Ref. 



III. THE PARTICLE VELOCITY FIELD 

In this section we discuss the compressibility of the 
particle velocity field. In particular, to clarify the range 
of validity of the Stokes approximation for the particle 
motion in fluid, in Sec. IIII Al we evaluate the particle re- 
sponse time T p as compared with the Kolmogorov micro- 
scale time t v . In Sec. IIII Bl we study the velocity of small 



and large particles in the turbulent fluid, required for 
evaluation of the effect of turbulent diffusion in these 
two regimes. This study also allows us to find the de- 
pendence of the degree of compressibility of the particle 
velocity field, o~ v , on the particle response time. 



A. Characteristic time scales and validly of the 
Stokes approximation 

Let us assume that the radius of the particles is small, 
so that the particle Reynolds number lZe p is smaller 
than the critical value, lZe cr , at which the laminar flow 
over a particle looses its stability. Then, we can ap- 
ply the Stokes approximation. It states that the fluid- 
particle friction force is proportional to the slip velocity, 
the difference between the particle velocity and the fluid 
velocity. A careful analysis performed by Lumley |2fij 
shows that in a turbulent flow the validity condition for 
lZe p < lZe cl - may be expressed via the particle radius a 
and the Kolmogorov micro-scale rj as follows: 



a < 2f 1 {p/ P p) 1 '\ 



(3.1) 



In the present analysis the ratio of the inertial time 
scale of the particles (the Stokes time scale t p ) and the 
turnover time of 77-eddies in the Kolmogorov micro-scale 
t(t)), is of primary importance. The particle response 
time is given by 



nip 2p p a 2 
p Qnvpa 9pv 

where the particle mass m p is: 

47T 

Y 



-a?p p 



(3.2) 



(3.3) 



The Kolmogorov micro-scale 77 is determined from the 
condition that the Reynolds number for eddies of scale r\ 
is equal to 1: 



lZe v = rjv(rj)/v = 1 



(3.4) 



Here v(rf) is the characteristic velocity of 77-scale eddies. 
It is related to the turnover time of the eddies as t(t)) = 
T}/v(rf). This allows us to rewrite the expression l|3.4l) as 
follows: 



r(n) = if /v. 



(3.5) 



Then, the ratio of the time-scales r p and r{rj) follows 
from Eqs. (22) and {53) : 



Tp _ 2ppO? 

r(rj) 9pn 2 



(3.6) 



Substituting the condition (|3.1I) for the validity of the 
Stokes approximation we find 



T p < ( Pp 



1/3 



(3.7) 
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Equation (|3.7|) implies that for "heavy" particles in a 
gas, that satisfy the Stokes approximation, the particle 
response time scale may be about ten times larger than 
the Kolmogorov time scale: t p < 10t(^). 



B. Particle velocity in turbulent fluid 

1. Equation for a particle velocity 

Assuming that the particles are small enough such that 
Eq. (|3.1|l is valid, we can apply Stokes' law for the fluid- 
particle friction force F p (t, r): 

F p (i,r)=C[v p (t,r)-u(t,r)], (3.8) 

with the Stokes friction coefficient C given by 

C = 6 7r pv a . (3-9) 

The Newton equation for a particle reads: 

dv p {t,r) 



dt 



-F p {t,r) = C[u(t,r)-v p (t,r)] . (3.10) 



Here the total time derivative (d/dt) takes into account 
the time dependence of the particle coordinate r: 



d 
dt 



d_ 

at 



v„(t,r). V 



(3.11) 



Now Eq. (|3.10(l takes the form: 



lU p (t,r) = u(i,r) . (3.12) 



|+v p (t,r).V 



In the following we analyze this equation in two limit- 
ing cases: for small particles with r p smaller than the 
turnover time of ^-eddies fSec. HII B 21) . and for large par- 
ticles, in Sec. IIII B 31 



2. Velocity of small particles 



Here the inverse operator is understood as its Taylor se- 
ries expansion: 

2 



duj (t,r) 
dt 



v e (t, r) = u4t,r) - t p — - — + 

(3.14) 

In the linear in t p approximation, this solution becomes 



vt(t,r) = ue(t,r) - t p 



u e {t,r) , 

(3.15) 

where in the RHS we replaced v e (t,r) — ► u e (t,r). Con- 
sequently, in this approximation 



div w(t,r) 



dr c 



at 



duf(t,r) du#(t, r) du"(t, r) 



dr? 



(3.16) 



(see Ref. ^3). Equation 1)3.16)) allows to determine the 
compressibility parameter a v , i)2.19)l . as follows: 



a v 



t(£) 



2 -t)' (;)'(?) 



4/3 



(3.17) 



In deriving this equation we estimated [du"(t, r)/dr a ] as 
[1/t(1)] and used Eqs. and l|2~5)l . 

Let a* be the characteristic size of particles for which 
t p — t(t)) for £ — r\ and, respectively, a v = 1: 

a*=T 1 yj9p/2pp . (3.18) 
Using this notation Eq. (|3.17|) can be rewritten as 

t p \ a J \r\) 

"(?r- <-») 



Equations 1)3.17(1 and 1)3.20(1 are valid only if o~ v < 1, 
otherwise the approximation of small t p /t(£) ratio is vi- 
olated. 



In this section we consider small particles for which 
the Stokes time t p is smaller than the turnover time of 
i'-eddies, t(£) sa £/u{£). These particles are completely 
involved in the motion of ^-eddies and for t p /t{£) = 0, we 
can write ve(t,r) = ug(t,r). Here subscript "£" denotes 
that we are dealing with the velocity field of ^-eddies, 
ue (t, r), and the velocity of particles in ^-clusters, vi(t, r), 
that was generated by f-eddies. In this approximation 
the velocity field ve(t, r) is incompressible. Therefore, 
the compressibility parameter a v may be determined by 
the first order corrections. 

To find the corrections to ~vi(t, r) linear in t p , we con- 
sider a formal solution of Eq. 1(3.12)1 in the form: 



Vf(£,r) 



d , 
T *dt + 1 



u e (t,r) 



(3.13) 



3. Velocity of large particles 

a. Effective equation of motion. In this section we 
consider the opposite case: the large particles with r p is 
larger than the turnover time of the smallest eddies in 
the Kolmogorov microscale t(j]), but smaller than the 
turnover time of the largest eddies t(L). Denote by £* 
the characteristic scale of eddies for which 







(3.21) 



This scale as well as the particles cluster size was intro- 
duced in Ref. 0. The eddies with £ ^> £„ almost fully 
involve particles in their motions, while the eddies with 
£ <§; £* do not affect the particle motions within the zero 
order approximation in the ratio [t(£)/t p ] <C 1. 
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To determine ve.it, r) we consider Eq. I|3.12|) in the co- 
moving with ^-eddies frame, where the surrounding fluid 
velocity u equals to the relative velocity of the ^-eddy at 
r, i.e., u(t,r) => Ui{t, r). At this point one has to take 
into account that the £-eAAy is swept out by all ^'-eddies 
of larger scales, £' > £, while the particle participates 
in motions of ^'-eddies with £' > £* > I. Therefore, 
the relative velocity Ug of the £-eddy and the particle 
is determined by ^'-eddies with the intermediate scales, 
£* > £' > £■ This velocity is determined by the con- 
tribution of £*-cddies, and can be considered as a time 
and space independent constant it* during the life time 
of the £ eddy and inside it. Velocity m» in our approach 
is random and we have to average the final result over 
statistics of ^*-eddies. Then Eq. (|3.12|l becomes 

d \ 

Tp dt + 1 j Vl<yt " V ' = r + U * t ' 

-r p [v e {t,r)- V\v t {t,r) . (3.22) 

In Eq. H3.22|l the velocity ug is calculated at point r and 
the velocity Vg is at r — n* t. For the sake of convenience 
we redefine here r — it* t = r' r and, respectively, 
r = r' + w» t => r + it* t. 

b. First non-vanishing contribution to vg. Clearly, 
vg(t,r) <C ug for £ <C £*■ Therefore we can find the 
first non- vanishing contribution to vg(t,r) in the limit 
[t(£) <C t p ] by considering the linear version of Eq. H3.22fl : 



d \ 

t p — + 1J v e (t, r) = ug(t, r + u*t) . (3.23) 



In the to, k representation this equation takes the form: 

{iujTp + 1) ve(u>, k) — ug(uj — k ■ u*,k) , (3.24) 

that allows one to find the relationship between the sec- 
ond order correlation functions F"g(ui, k) and F"^(ui, k) 
of the velocity fields vg and ug: 



1 F°f(w-fc-u„fc) . (3.25) 



w 2 r p 2 + l u 



Functions F"%(u, k) and F^(u>, k) are defined as usual: 

(2^) 4 <5(w + u>') S(k + k')FZ$(u, k) (3.26) 
= (vf(w,k)vjl(u',k'jS , etc. 

The simultaneous correlation functions are related to 
their w-dependent counterparts via the integral / dw/27r, 
e.g., 



(3.27) 



The tensorial structure of F^(k) follows from the incom- 
pressibility condition and the assumption of isotropy: 



F^(k)=P^(k)F u4 (k), 



(3.28) 



where P Q/3 (fc) is the transversal projector: 
P a ^(k) = S al3 -k a k /k 2 . 



(3.29) 



In the inertial range of scales the function F"^(w, k) may 
be written in the following form: 

F^(uj, k) = P a P{k) F u j(k) t(£) f[u> t(£)} . (3.30) 

Here the dimensionless function f(x) is normalized as 
follows: 



fix) dx — 2-7T 



(3.31) 



Now we can average Eq. (|3.25|l over the statistics of £*- 
eddies. Denoting the mean value of some function gix) 
as gix) we have: 



/[(w-fe-u,)T(/)] » 5[(w - k ■ it*) t(£)] (3.32) 

£ fu*<jj 
J- 



t{1) u* 



£ 



Here the dimensionless function /* (x) has one maximum 
at x = 0, and it is normalized according to Eq. I|3.31|l . 
The particular form of /*(x) depends on the statistics of 
^*-eddies and our qualitative analysis is not sensitive to 
this form. Thus, we may choose, for instance: 



Uix) = 2/[x 2 + 1] . 



(3.33) 



In Eqs. (|3.32|) we took into account that the characteristic 
Doppler frequency of £-eddies (in the random velocity 
field m* of £*-eddies) may be evaluated as: 



(3.34) 



This frequency is much larger than the characteristic fre- 
quency width of the function /[wt(Q] (equal to l/r(£)), 
and therefore the function fix) in Eq. (|3.32|) may be ap- 
proximated by the delta function Six). 

After averaging, Eq. 13.25f) may be written as 



P a fl(fc)/.(0) £ 

UJ 2 T 2 + 1 U* 



F u 4k) 



(3.35) 



Here we took into account that t p 3> £/u* that allows us 
to neglect the frequency dependence of /*(u* u>/£) and to 
calculate this function at u> = 0. Together with Eq. (|3.35|) 
this yields 



Ffgik) = P af \k)F u Ak) 



a/3. 



(3.36) 



where we used the estimate /*(0) ~ 2, that follows from 

Eq. ESSl. 

The equation (|3.36(l provides the relationship between 
the mean square relative velocity of ^-separated particles, 
vg, and the velocity of f -eddies, ug: 



£ £ 

Vg ~ Ug, ~ Ug\ — 

y Tp-u, v £* 



(3.37) 
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c. Effective nonlinear equation. For a qualitative 
analysis of the role of the nonlinearity of the particle 
behavior in an ^-cluster we evaluate V in the nonlinear 
term, Eq. ((3.221) . as l/£, neglecting the spatial depen- 
dence and the vector structure. The resulting equation 
in cj-representation reads: 



(ioj + ') p )Vi(u)=^ Al U l (u)+]\f u) , 7p = V' 



p ■ 



Vi(uj) — J ve(t) exp[— iujt]dt , 
1 f 

v e(t) — tt I V^(uj) exp[iu}t]dui . 



(3.38) 



In the zeroth order (linear) approximation [M w — > 0) 

(3.39) 



which is the simplified version of Eq. (|3.24|) . This allows 
us to find in the linear approximation 



Mt)) = J = J ^-t—t , (3-40) 

where F u j(uj) is the correlation functions of Ue(ui): 

2wS(uj + u')F u ,t{u) = (Ut(u)U t (u')) , (3.41) 

similarly to Eq. 1)3.26(1 . 

In the limit t p ^> £/u* one can neglect in Eq. I|3.4U|) the 
w-dependence of F u ^{ui), which has characteristic width 
l/ust and conclude: 



vi fi ee (vi(t))^^-F u A0)^uf-^^ui- : 



.2 _ /„,2 



uf(t)) 



(3.42) 



in agreement with Eq. (|3.37l) . 

d. First nonlinear correction. To evaluate the first 
nonlinear correction to Eq. I|3.42|) one has to substitute 
V £ (uj) from Eq. $nty into Eq. (|3~33|) for Af u : 



7 

V^x{uj) = I du!idLU 2 S(uj + tvi + oj 2 ) 

Ue(uji) Ui(lo 2 ) 



lUi + 7 P iuj 2 + 7 p 



(3.43) 



Using Eq. (|3.43|l instead of Eq. I|3.39|l we obtain instead 
of Eq. (EHIl 

du 
M' 

diO\dw 2 



F u js(oj) 



= (N,iW] 2 > - / ^F uA1 {u), (8.441 



(o> 2 + 7 2)£2 



2tt 



F u j.{ui) F Ut i(u} 2 )6(u) + mi + uj 2 ) 



In this derivation we assumed for simplicity the Gaussian 
statistics of the velocity field. 
Now let us estimate 



( 2 



2/3 



(3.45) 



that is much larger than the result (|3.42|l for i; 2 ob- 
tained in the linear approximation. This means that the 
simple iteration procedure we used is inconsistent, since 
it involves expansion in large parameter [(£*/ 'i) 1 ^ 3 }. 

e. Renormalized perturbative expansion. A similar 
situation with a perturbative expansion occurs in the 
theory of hydrodynamic turbulence, where a simple it- 
eration of the nonlinear term with respect to the lin- 
ear (viscous) term, yields the power series expansion in 
IZe 2 3> 1. A way out, used in the theory of hydro- 
dynamic turbulence is the Dyson- Wyld re-summation 
of one-eddy irreducible diagrams (for details see, e.g., 
Refs. |2(| 153, HU). This procedure corresponds to ac- 
counting for the nonlinear (so-called "turbulent" viscos- 
ity) instead of the linear, kinematic viscosity. A similar 
approach in our problem implies that we have to account 
for the self-consistent, nonlinear renormalization of the 
particle frequency j p T p {£) in Eq. I|3.38|l and to sub- 
tract the corresponding terms from Af u . With these cor- 
rections, Eq. (|3.38|l reads: 



[iu + r p (£)]Vt(u>)=ipUt(u)+tf u 



(3.46) 



Here M u is the nonlinear term J\f u after substraction of 
the nonlinear contribution to the difference 



r P W-7p 



(3.47) 



The latter relation actually follows from a more detailed 
perturbation diagrammatic approach. In our context it 
is sufficient to realize that in the limit T p (£) 3> "f p one 
may evaluate T p (£) by a simple dimensional reasoning: 



(3.48) 



which is consistent with Eq. (|3.47|) . In addition, 
Eq. 13.47(1 has a natural limiting case T p (£) — > 7 P when 
vi/£ < 7 P . Now using Eq. l|3~l6} instead of Eq. jSHSj 
we arrive at: 



iLu + T p (£) ■ 

Accordingly, instead of the estimates 13.42(1 one has 



(3.49) 



J e,o 



,2 7 P 



T p (£) it* 



T p {£) \£ 



(3.50) 



The latter equation together with Eq. ((3.48(1 allows to 
evaluate T p (£) as follows: 



T P (£) 



2 2 

*4n 



1/3 



7 P 



1/9 



(3.51) 
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Hence the estimate (|3.50f) becomes 





\r(ty 


5/3 




. T P _ 





(3.52) 



Repeating the evaluation of the nonlinear correction v\ 2 
with the renormalized Eq. I|3.46fl we find that 



via 



(3.53) 



This means that now the expansion parameter is of the 
order of 1, in accordance with the renormalized pertur- 
bation approach. 



IV. THE CLUSTERING INSTABILITY OF THE 
SECOND MOMENT OF PARTICLE NUMBER 
DENSITY 

In this section we will present a quantitative analysis 
for the clustering instability of the second moment of 
particles number density. 



The form of the coefficients B(R), U(R) and D a p(R) 
depends on the model of turbulent velocity field. For 
instance, for the random velocity with Gaussian statistics 
of the Lagrangian trajectories ^(r |t) these coefficients are 
given by 

/>oo 

B(R) » 2 / (6[0,C(n|0)]6[T,€(r a |r)])dT, (4.3) 
Jo 

/>oo 

U(R) « -2 / («[0,€(ri|0)]6[r,€(r2|r)])dr, 
Jo 

/>oo 

D T af3 (R) » 2 / (v a [0,£{ri\0)]v fi [T,e(r*\T)])dT 
Jo 

(for details see Ref. where b = divv. 

B. Clustering instability 

Let us study the clustering instability. Consider the 
range of scales a < £ <C (*, where the size of a particle 
a>rj. Then the relationship between vf and uj reads: 



A. Basic equations 

To determine the growth rate of the clustering insta- 
bility let us consider the equation for the two-point cor- 
relation function $(<, R) of particle number density: 



<9$ 
~dt 



[B(R) + 2U(R) • V 



D af ,(R)V a V{,]*(t,R), 

(4.1) 

(see Ref. 0). The meaning of the coefficients B(R), 
U(R) and D a p(R) is as follows: 

The function B(R) is determined by the compressibil- 
ity of the velocity field and it causes the generation of 
fluctuations of the number density of particles. 

The vector U(R) determines a scale-dependent drift 
velocity which describes a transport of fluctuations of 
particle number density from smaller scales to larger 
scales, i.e., in the regions with larger turbulent diffusion. 
The latter can decrease the growth rate of the clustering 
instability. Note that U(i? = 0) = whereas B(R = 
0) 7^ 0. For incompressible velocity field V(R) = and 
B{R) = 0. 

The scale-dependent tensor of turbulent diffusion 
D a p(R) is also affected by the compressibility. In very 
small scales this tensor is equal to the tensor of the molec- 
ular (Brownian) diffusion, while in the vicinity of the 
maximum scale of turbulent motions this tensor coincides 
with the regular tensor of turbulent diffusion. 

Thus, the clustering instability is determined by the 
competition between these three processes. The tensor 
D a p(R) may be written as 



D a p(R) = 2DS aP + DJR 



■.p\ 

- T 



(4.2) 



D af3 (R) = D af} (0) - D af} (R) 



u 2 



t.o 



r(£) 



(4.4) 



where according to Eq. (|3.52|) the exponent s = 5/3. In 
this case the expression for the turbulent diffusion tensor 
in nondimensional form reads 

D T a JR) = R^^/^dRH^ + C 2 R a R p ) , (4.5) 



C x = 



C 2 - 



5 + 4s + 6cr T 

9(1 + ct t ) ' 
(4 8 -l)(2<7 T - 
9(l + cr T ) 



where R is measured in the units of L and time t is mea- 
sured in the units of t l = t(1 — L). The parameter a T 
is defined by analogy with Eq. 12.19fl : 



V D T V 

V x D„ x V 



V a V p DJR) 



(4.6) 



where e Q( 3 7 is the fully antisymmetric unit tensor. Equa- 
tions H2.19fl and (|4.6I) imply that <r T = a v in the case of 
(^-correlated in time compressible velocity field. 

For a random incompressible velocity field with a fi- 
nite correlation time the tensor of turbulent diffusion 
D T a p(R) = r- 1 (^(ri|i)^(r 2 |t)) (see Ref. E3) and the 
degree of compressibility of this tensor is 



((v-O 2 ) 

<(Vx£) 2 ) 



(4.7) 



where £(r|t) is the Lagrangian trajectory. 

To determine the functions B(R) and U(R) we use 
the general form of the two-point correlation function 
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of the particle velocity field in the the range of scales 

(v a (t,r)vp(t + r,r + R)} = ^[5 a p - {C v l R 2 5 a0 

+ C v 2 R a R p )R 2{ - s -^^]f{r), (4.8) 
(4 + 5 + 3^) 
1 3(l + a„) ' 

(1 + S )(2<7„ - f) 

3(l + o-„) 

Substitution Eq. I|4.8[) into Eq. (|4.3[) yields 

U(R) = U Q R^ 7)/3 , 
B{R)=B i?( 4s - ? )/ 3 , 

where 



ci = 



(4.9) 



f 



B = aU 



and the coefficients a and /3 depend on the properties of 
turbulent velocity field. The dimensionless functions B 
and Uo in Eq. 1)4. 9fl are measured in the units of r^ 1 . 

For the (5-correlated in time random Gaussian com- 
pressible velocity field a T — a v and 



B{R)=V a VpD aP {R). 
U a (R) = VpD af3 (R) 



(4.10) 
(4.11) 



(for details see Ref. dill El). In this case the second 
moment $(i, R) can only decay, in spite of the compress- 
ibility of the velocity field. 

For the finite correlation time of the turbulent velocity 
field <7 T a v and Eqs. I|4.10|l and (|4.11|) are not valid. 
The clustering instability depends on the ratio a T /a v . 
In order to provide the correct asymptotic behaviour of 
Eq. H4.9fl in the limiting case of the 5-correlated in time 
random Gaussian compressible velocity field we have to 
choose the coefficients (3 and a in the form: 

/? = 8(4s 2 + 7s-2)/27, a=(4s + 2)/3. 

Note that when s < 1/4 the parameters (3 < and 
B(R) < 0. In this case there is no clustering instabil- 
ity of the second moment of particle number density. 
Thus, Eq. H4.1|) in a non-dimensional form reads: 



^ = i?( 4s - 7 )/ 3 [i? 2 $"(d 
dt 

+B <S>] . 



C 2 ) + 2R&(U + Cx) 

(4.12) 



Consider a solution of Eq. H4.12|l in the vicinity of the 
thresholds of the excitation of the clustering instability, 
wh ere (d <^/dt)R^ 7 - 4s ^ 3 is very small. Thus, the solution 
of l|4.12|) in this region is 



$(i?) = A X BT 



(4.13) 



where Ai : 
A 



A ± 

Ci — C 2 



2U Q 



C 3 



V = 



2(Cx+C 2 ) ' " 2(C X +C 2 )' 
C\ = 45 (C 1 + C 2 ) - (d -C 2 + 2U ) 2 . 

Since the correlation function Q(R) has a global maxi- 
mum at R = a, the coefficient C\ > C 2 — 2[/ if // is 
a real number (see below). Thus the asymptotic solu- 
tion of the equation for the two-point correlation func- 
tion <f>(t, R) of the particle number density in the range 
of scales a < £ < I* reads 



where 



$(i?) = Aiir A sin[/zlii(4/#)] > 



(4.14) 



a ) sin[/xln(^*/a)] 



Now consider the range of scales f, <C ^ < £■ In this 
case the nondimensional form of the turbulent diffusion 
tensor is given by 

D ap {R) = i?- 2 / 3 (Cii? 2 ^ + C 2 R a R p ) , (4.15) 
2(5 + 3a T ) - 4(2<r T -l) 



C x = 



9(1 



C 2 



9(1 



and Eq. (|4.1|1 reads: 



<9$ 
~dt 



R- 2 ^[R 2 ^"{Cx + C 2 ) +2R<f>'Cx] . (4.16) 



Here we took into account that in this range of scales the 
functions B(R) and U{R) are negligibly small, and the 
degree of compressibility of the turbulent diffusion tensor 
can be different in the above two ranges of scales. 

Consider a solution of Eq. (|4.1(j|) in the vicinity of 
the thresholds of the excitation of the clustering instabil- 
ity, when (d<& / dt)R 2 / 3 is very small. Thus, the solution 
of (|4.16(l is given by 



where 



$(R) = A 2 R- 



\Cx - C 2 



(4.17) 



Cx + C 2 



7cL 



The growth rate of the second moment of particle num- 
ber density can be obtained by matching the correlation 
function & (R) and its first derivative $'(i?) at the bound- 
ary of the above two ranges of scales, i.e., at the points 
R — 1*1 L. The matching yields that in the range of scales 

-C I < L the asymptotic solution of the equation for 
the two-point correlation function <I>(i, R) has the form 
of Eq. ijCTfjl with 

a 2 -a+i 

A 2 = (-l) k Ax(^ 
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FIG. 1: The range of parameters (a v ,a T ) for which the 
clustering instability may occur. The various curves indi- 
cate results for s = 7/4 (dashed), s — 5/3 (solid), for s = 1 
(dashed-dotted) and for s = 2/3 (dotted). The thin dashed 
line a v = <r T corresponds to the ((-correlated in time random 
compressible velocity field. 



Such matching is possible only when Ai is a complex 
number, i.e., when Cf > (i.e., /i is a real number). The 
latter determines the necessary condition for the cluster- 
ing instability of particle spatial distribution. The range 
of parameters (a v , a T ) for which the clustering instabil- 
ity of the second moment of particle number density may 
occur is shown in Fig. 1. The line a v = er T corresponds 
to the ^-correlated in time random compressible veloc- 
ity field for which the clustering instability cannot be 
excited. The various curves indicate results for different 
value of the parameter s. In particular, the value s = 7/4 
corresponds to the turbulent diffusion tensor with the 
scaling cx R 2 [see Eq. 1)4. 5(1 ]. The curves for s = 7/4 
(dashed) and s = 5/3 (solid) practically coincide. The 
parameter s can be considered as a phenomenological pa- 
rameter, and the change of this parameter from s = 5/3 
to s = can describe a transition from one asymptotic 
behaviour (in the range of scales a < £ < ) to the other 
(4 <•£<£)• 



V. DISCUSSION 

Formation and evolution of particle clusters are of 
fundamental significance in many areas of environmen- 
tal sciences, physics of the atmosphere and meteorol- 
ogy (e.g., smo g a nd fog formation, rain formation), see 
e.g., Ref. [13, |33, El H3, transport and mixing in in- 
dustrial turbulent flows (like spray drying and cyclone 
dust separation, dynamics of fuel droplets), see e.g., 



Ref. |33j |3J, |35J, |3jJ |3J, |3J and references therein. The 
analysis of the experimental data showed that the spatial 
distributions of droplets in clouds are strongly inhomo- 
geneous |l2j . The small-scale inhomogeneities in particle 
distribution were observed also in laboratory turbulent 

flows E3G3- 

The analyzed effect of particle clustering may be of 
relevance in turbulent fluid flows of different nature with 
inertial particles or droplets (e.g., in atmospheric turbu- 
lence, combustion and in a laboratory turbulence). In 
particular, this effect can cause formation of small-scale 
inhomogeneities in spatial distribution of fuel droplets in 
diesel engines. The starting point is the above theoreti- 
cal model that describes the inertial particle clustering. 
For numerical estimates we adopt the operating param- 
eters for a typical diesel engines, taking the crankshaft 
rotational speed iV=3000 rev/min, that corresponds to 
the average piston speed, v p = 10 3 cm/s, and the turbu- 
lent velocity v' ~ 5 • 10 2 cm/s. Taking into account that 
the spatial integral space (the characteristic scale of the 
largest turbulent eddies) for this conditions is L ~ 0.3 
cm (see Ref. [1^), we find for the Reynolds number: 



Tie 



v'L 



10 4 



with 



(5.1) 



5 • 10 2 cm/sec, L~ 0.3 cm 



Here the kinematic viscosity of fluid is v ps 0.015 cm 2 /s. 
To estimate v we used the following values for the fluid 
density p = 0.025 g/cm 3 (at 62 atm. and T=860 K), and 
for kinematic viscosity 



fl(3-7)/2 



(5.2) 



where R = (p/po) ~ 15 is the compression ratio. Then 
the Kolmogorov micro-scale is r\ « L/TZe 3 ^ 4 ~ 3 pm. A 
more accurate estimate is 



n 



{Ke/Ke^fl 4 ' 



(5.3) 



where TZe cr > 1 is the critical value of TZe at which the 
laminar flow becomes unstable. Taking lZe CT i=s 5, we 
obtain a more realistic value rj ~ 10 pm. 

Note that turbulence in the in diesel engines is neither 
homogeneous nor isotropic. However, we study dynam- 
ics of droplets at very small scales which are of the order 
of tens of hundreds microns. In such small scales the 
turbulence can be considered as a quasi isotropic and 
homogeneous. Anisotropy and inhomogeneities of turbu- 
lence can be essential in scales which are of the order of 
1 cm and larger. 

Consider a region of unstable cluster scales: rj < 
£ < £*. Assuming the fuel density in droplets pdr = 
0.85 g/cm 3 , we find from the expression l|3.18|l the esti- 
mate for the critical value of the droplet radius, a* ~ 
4 pm, indicating that all droplets with a > 4 pm are un- 
stable against formation of the clusters with scales rang- 
ing from rj to This equation yields the following val- 
ues of £*: is, ~ 160 pm for a = 10 pm and ~ 1cm for 
a — 40 pm. 
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The characteristic time of the cluster growth t c i(£) — 
1/jd depends on the cluster size i. Estimate for £ = £„/3 
gives r c i(^*/3) « 1.7t{£*) Td r . The turnover time t{£) 
can be estimated as follows: 



r(£) 



2/3 , g x 2/3 

~3xlO~ 3 (-J sec. (5.4) 



Here we used our estimates l|5.1|l for the integral scale L 
and turbulent velocity v' . Assuming for the droplets with 
the size a = 10 /im, £ = £*/3 w 50 pm (with £„ = 160 pm) 
one has t c \{£) ~ 0.35 /us. For the droplets with the radius 
a = 40 /im with £* w 1 cm, one has r c i(£*/3) ~ 5.5 ps, 
which is larger than the characteristic spray time in the 
combustion chamber (a few milliseconds). These esti- 
mates imply that for the chosen parameters the distri- 
bution of the droplets in the spray will be shifted to the 
droplets with larger diameters at the given turbulent in- 
tensity. 

Now we estimate the time of the turbulent air- fuel mix- 
ing, r T ~ 1/t t in the combustion chamber. Taking into 
account that 



7t 



~~dJ 



v'L 
3d? 



(5.5) 



where D T is the coefficient of turbulent diffusion, and 
d ~ 10 cm is the bore diameter, we find for the time of 
turbulent air-fuel mixing r T » Is. Thus, the turbulent 
diffusion itself is too slow and a strong external global 
flow - "swirl" is needed for the effective air-fuel mixing. 

In the previous sections we have shown that the parti- 
cle spatial distribution in the turbulent flow field is un- 
stable against formation of clusters with particle num- 
ber density that is much higher than the average parti- 
cle number density. Obviously this exponential growth 
at the linear stage of instability should be saturated by 
nonlinear effects. The nonlinear saturation for the water 
droplets in the atmospheric turbulence (in clouds) was 
discussed in Ref . 01 • Here we discuss a similar issue for 
the fuel droplets under the typical conditions pertinent 
to diesel engines. 

A momentum coupling of particles and turbulent fluid 
is essential when m p n c \ ~ p, i.e., the mass loading pa- 
rameter, <j) = m p n c \/ p, is of the order of unity (see, e.g., 
Ref. H||). This condition implies that the kinetic energy 
of air p(u 2 ) is of the order of the droplets kinetic energy 
w p n c i(v 2 ), where |u| ~ |v|. This yields: 



n c \ ~ a (p/3p p ) 



(5.6) 



For the fuel droplets in the diesel combustion chambers 
(see, e.g. Ref. yij) p P /p — 34. Thus, e.g., for a ~ 10 pm 
we obtain n c \ ~ 3 x 10 7 cm -3 . For the cluster with the 
size £ w £*/3 « 50 pm this yields for the total number of 
particles in the cluster of that size N c \ ~ £ 3 n c \ ~ 30. Note 
that the mean number density of droplets in a combus- 
tion camera n is about 10 4 cm -3 . Therefore the clustering 



instability of droplets in the diesel engines increases their 
concentrations in the clusters by the orders of magnitude. 

Note that in the initial stage of the clustering instabil- 
ity, when the mass loading parameter is small, i.e., for 
the small number density of droplets, we can neglect the 
droplets collisions, and consider only the collisions be- 
tween the droplets and the molecules of the surrounding 
fluid (air). In the nonlinear stage of the clustering in- 
stability, the four-way coupling can be effective, and the 
kinetic approach for analysis of the droplets collisions 
inside the cluster can be important. Moreover, the co- 
agulations of droplets due to their collisions change the 
size of the droplets. This effect is determined by the 
kinetic Smoluchovsky equation for droplets size distri- 
bution. However, the interaction between droplets and 
the surrounding fluid even at this stage of the clustering 
instability can be described on the level of the contin- 
uum hydrodynamic approach. The main purpose of the 
present paper is to describe the initial (linear) stage of 
the clustering instability and to determine the conditions 
for the onset of the clustering instability of fuel droplets. 
This stage of the instability can be analyzed only us- 
ing the continuum hydrodynamic approach because the 
mean free path l c of the molecules of the surrounding fluid 
is much smaller than all the characteristic spatial scales in 
the problem (e.g., the Kolmogorov micro-scale 77, the size 
of droplets a, etc). Indeed, for diesel engines these param- 
eters are: l c = 0.01 pm, 77 = 3 pm, a > 10 pm. Similarly, 
the mean time r c between collisions of the molecules of 
the surrounding fluid is much smaller than all the char- 
acteristic time scales (e.g., the correlation time in the 
Kolmogorov micro-scale t v , the Stokes time for droplets 
r p , etc). Indeed, for diesel engines these parameters are 
t c = 2 x 10- n s, t v = 3 x 10~ 5 s, and r p = 10~ 4 s for 
a = 10 pm. At the strongly nonlinear stage of the cluster- 
ing instability, when the mean free path of the droplets 
inside a cluster becomes much smaller than the size of 
the cluster, the continuum hydrodynamic approach can 
be still used to obtain a rough estimate of the number 
density of droplets inside the cluster at the saturation of 
the clustering instability. 
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